Packages
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") # workflow and plots
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # for p-values
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("UsingR")) install.packages("UsingR")
library("UsingR") # simple.eda model assumption checker
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # lmer model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans') # marginal means
if (!require("AICcmodavg")) install.packages("AICcmodavg")
library("AICcmodavg") # model selection
if (!require("MuMIn")) install.packages("MuMIn")
library("MuMIn") # model selection
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # color
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown formatBackground and Goals
This data was collected June - August by Master’s student Savannah Weaver, advisor Dr. Emily Taylor, and research assistants Tess McIntyre and Taylor Van Rossum. Adult male Sceloporus occidentalis were caught across the Cal Poly campus then acclimated to 4 different climate treatments. This R file analyzes the effect of experimental climate treatments on lizard body condition, osmotic balance (plasma osmolality and hematocrit), and cutaneous evaporative water loss (CEWL).
Data
Load
Read-in data that was compiled, formatted, and checked for completeness in ‘wrangling_general’. See that file for information related to the variables.
## measurement_date type individual_ID mass_g
## Min. :2021-06-16 exp :804 201 : 7 Min. : 7.00
## 1st Qu.:2021-07-01 rehab:132 202 : 7 1st Qu.: 9.50
## Median :2021-07-25 203 : 7 Median :10.60
## Mean :2021-07-22 204 : 7 Mean :10.64
## 3rd Qu.:2021-08-14 205 : 7 3rd Qu.:11.70
## Max. :2021-09-01 206 : 7 Max. :17.40
## (Other):894
## hematocrit_percent trial_number temp_tmt humidity_tmt SVL_mm
## Min. :13.00 1:175 Hot :467 Humid:468 Min. :60.00
## 1st Qu.:26.00 2:203 Cool:469 Dry :468 1st Qu.:66.00
## Median :32.00 3:231 Median :67.00
## Mean :31.99 4:189 Mean :67.74
## 3rd Qu.:37.00 5:138 3rd Qu.:70.00
## Max. :52.00 Max. :77.00
## NA's :408
## tmt day_n day_factor osmolality_mmol_kg_mean
## Cool Humid:238 Min. : 0.000 0 :134 Min. :295.3
## Hot Humid :230 1st Qu.: 4.000 4 :134 1st Qu.:336.1
## Cool Dry :231 Median : 6.000 5 :134 Median :351.3
## Hot Dry :237 Mean : 5.705 6 :134 Mean :354.3
## 3rd Qu.: 8.000 7 :134 3rd Qu.:370.0
## Max. :10.000 8 :134 Max. :471.5
## 10:132 NA's :414
## CEWL_g_m2h_mean msmt_temp_C msmt_RH_percent cloacal_temp_C
## Min. : 7.152 Min. :24.80 Min. :25.52 Min. :23.00
## 1st Qu.:19.755 1st Qu.:26.20 1st Qu.:46.11 1st Qu.:25.00
## Median :24.152 Median :26.74 Median :47.88 Median :26.00
## Mean :24.767 Mean :26.72 Mean :46.74 Mean :25.92
## 3rd Qu.:28.505 3rd Qu.:27.11 3rd Qu.:50.50 3rd Qu.:27.00
## Max. :56.066 Max. :29.20 Max. :56.16 Max. :30.00
## NA's :669 NA's :668 NA's :668 NA's :668
## msmt_temp_K e_s_kPa_m e_a_kPa_m msmt_VPD_kPa
## Min. :297.9 Min. :3.219 Min. :0.9894 Min. :1.486
## 1st Qu.:299.4 1st Qu.:3.504 1st Qu.:1.6464 1st Qu.:1.767
## Median :299.9 Median :3.620 Median :1.7411 Median :1.853
## Mean :299.9 Mean :3.620 Mean :1.6833 Mean :1.937
## 3rd Qu.:300.3 3rd Qu.:3.701 3rd Qu.:1.7992 3rd Qu.:2.012
## Max. :302.4 Max. :4.194 Max. :1.9326 Max. :3.021
## NA's :668 NA's :668 NA's :668 NA's :668
## SMI temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
## Min. : 6.747 Min. :23.30 Min. :0.5966 Min. :13.75
## 1st Qu.: 9.714 1st Qu.:24.05 1st Qu.:0.7828 1st Qu.:29.21
## Median :10.594 Median :24.88 Median :1.0461 Median :45.24
## Mean :10.599 Mean :29.60 Mean :1.1513 Mean :52.94
## 3rd Qu.:11.390 3rd Qu.:35.05 3rd Qu.:1.5191 3rd Qu.:82.84
## Max. :15.063 Max. :36.00 Max. :1.8447 Max. :93.15
##
## humidity_SD_tmttrial e_a_mean_tmttrial e_a_SD_tmttrial VPD_mean_tmttrial
## Min. : 4.370 Min. :0.4258 Min. :0.1891 Min. :0.2053
## 1st Qu.: 6.234 1st Qu.:1.3896 1st Qu.:0.2651 1st Qu.:0.8018
## Median : 7.382 Median :2.2917 Median :0.3166 Median :2.0450
## Mean : 8.765 Mean :2.3340 Mean :0.3936 Mean :2.0043
## 3rd Qu.:12.297 3rd Qu.:2.6749 3rd Qu.:0.5019 3rd Qu.:3.1651
## Max. :19.846 Max. :4.8058 Max. :0.7691 Max. :4.0685
##
## VPD_SD_tmttrial temp_C humidity_percent e_a_kPa
## Min. :0.1690 Min. :23.80 Min. :17.90 Min. :0.50
## 1st Qu.:0.2558 1st Qu.:23.80 1st Qu.:34.10 1st Qu.:2.00
## Median :0.3483 Median :24.40 Median :56.20 Median :2.15
## Mean :0.4038 Mean :29.64 Mean :52.84 Mean :2.32
## 3rd Qu.:0.4735 3rd Qu.:35.50 3rd Qu.:78.30 3rd Qu.:2.30
## Max. :0.8863 Max. :35.50 Max. :80.90 Max. :4.50
##
## VPD_kPa
## Min. :0.600
## 1st Qu.:0.600
## Median :1.800
## Mean :2.002
## 3rd Qu.:3.800
## Max. :3.800
## Split
Make sub-dataframes without recovery data / with only recovery-related data:
note: recovery analysis was not included in publication
dat_no_rehab <- dat %>%
dplyr::filter(day_n %in% c(seq(0,8)))
recovery_values <- dat %>%
dplyr::filter(day_n == 10) %>%
dplyr::select(individual_ID,
end_hct = hematocrit_percent,
end_osml = osmolality_mmol_kg_mean,
end_SMI = SMI)
recovery_v_post_exp <- dat %>%
dplyr::filter(day_n == 8) %>%
left_join(recovery_values, by = 'individual_ID') %>%
mutate(delta_osml_10_8 = end_osml - osmolality_mmol_kg_mean,
delta_hct_10_8 = end_hct - hematocrit_percent,
delta_SMI_10_8 = end_SMI - SMI)
recovery_v_pre_exp <- dat %>%
dplyr::filter(day_n == 0) %>%
left_join(recovery_values, by = 'individual_ID') %>%
mutate(delta_osml_10_0 = end_osml - osmolality_mmol_kg_mean,
delta_hct_10_0 = end_hct - hematocrit_percent,
delta_SMI_10_0 = end_SMI - SMI)Check
Dates:
## [1] "2021-06-16" "2021-06-20" "2021-06-21" "2021-06-22" "2021-06-23"
## [6] "2021-06-24" "2021-06-26" "2021-06-30" "2021-07-01" "2021-07-02"
## [11] "2021-07-03" "2021-07-04" "2021-07-06" "2021-07-20" "2021-07-24"
## [16] "2021-07-25" "2021-07-26" "2021-07-27" "2021-07-28" "2021-07-30"
## [21] "2021-08-08" "2021-08-12" "2021-08-13" "2021-08-14" "2021-08-15"
## [26] "2021-08-16" "2021-08-18" "2021-08-22" "2021-08-26" "2021-08-27"
## [31] "2021-08-28" "2021-08-29" "2021-08-30" "2021-09-01"Number of measurements for each lizard:
## # A tibble: 134 × 2
## individual_ID n
## <fct> <int>
## 1 201 6
## 2 202 6
## 3 203 6
## 4 204 6
## 5 205 6
## 6 206 6
## 7 207 6
## 8 208 6
## 9 209 6
## 10 210 6
## # … with 124 more rowsEvery lizard has 6 experimental measurements: pre-tmt, mid-tmt, post-tmt, and mass checks on each of the 3 days between mid and post-tmt.
Did any of the treatment groups inherently start out with large differences in response variables?
dat %>%
dplyr::filter(day_n == 0) %>%
group_by(tmt) %>%
summarise(mean(mass_g),
sd(mass_g),
mean(SMI),
mean(hematocrit_percent),
mean(osmolality_mmol_kg_mean),
mean(CEWL_g_m2h_mean))## # A tibble: 4 × 7
## tmt `mean(mass_g)` `sd(mass_g)` `mean(SMI)` mean(hema…¹ mean(…² mean(…³
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cool Humid 11.6 1.35 11.7 39.6 351. 20.9
## 2 Hot Humid 11.6 1.75 11.5 37.9 347. 21.4
## 3 Cool Dry 11.8 1.61 11.8 39.3 346. 20.0
## 4 Hot Dry 12.0 1.68 11.8 38.9 347. 20.9
## # … with abbreviated variable names ¹`mean(hematocrit_percent)`,
## # ²`mean(osmolality_mmol_kg_mean)`, ³`mean(CEWL_g_m2h_mean)`There are slight differences, but overall the starting values across groups are more or less the same.
Temp & RH during (all, before and after exp) CEWL measurements:
## measurement_date type individual_ID mass_g
## Min. :2021-06-16 exp :804 201 : 6 Min. : 7.00
## 1st Qu.:2021-06-30 rehab: 0 202 : 6 1st Qu.: 9.50
## Median :2021-07-25 203 : 6 Median :10.60
## Mean :2021-07-22 204 : 6 Mean :10.65
## 3rd Qu.:2021-08-13 205 : 6 3rd Qu.:11.60
## Max. :2021-08-30 206 : 6 Max. :17.40
## (Other):768
## hematocrit_percent trial_number temp_tmt humidity_tmt SVL_mm
## Min. :13.00 1:150 Hot :402 Humid:402 Min. :60.00
## 1st Qu.:28.25 2:174 Cool:402 Dry :402 1st Qu.:66.00
## Median :33.00 3:198 Median :67.00
## Mean :33.75 4:162 Mean :67.73
## 3rd Qu.:39.00 5:120 3rd Qu.:70.00
## Max. :52.00 Max. :77.00
## NA's :406
## tmt day_n day_factor osmolality_mmol_kg_mean
## Cool Humid:204 Min. :0.0 0 :134 Min. :295.3
## Hot Humid :198 1st Qu.:4.0 4 :134 1st Qu.:334.7
## Cool Dry :198 Median :5.5 5 :134 Median :348.3
## Hot Dry :204 Mean :5.0 6 :134 Mean :352.3
## 3rd Qu.:7.0 7 :134 3rd Qu.:367.4
## Max. :8.0 8 :134 Max. :445.5
## 10: 0 NA's :413
## CEWL_g_m2h_mean msmt_temp_C msmt_RH_percent cloacal_temp_C
## Min. : 7.152 Min. :24.80 Min. :25.52 Min. :23.00
## 1st Qu.:19.755 1st Qu.:26.20 1st Qu.:46.11 1st Qu.:25.00
## Median :24.152 Median :26.74 Median :47.88 Median :26.00
## Mean :24.767 Mean :26.72 Mean :46.74 Mean :25.92
## 3rd Qu.:28.505 3rd Qu.:27.11 3rd Qu.:50.50 3rd Qu.:27.00
## Max. :56.066 Max. :29.20 Max. :56.16 Max. :30.00
## NA's :537 NA's :536 NA's :536 NA's :536
## msmt_temp_K e_s_kPa_m e_a_kPa_m msmt_VPD_kPa
## Min. :297.9 Min. :3.219 Min. :0.9894 Min. :1.486
## 1st Qu.:299.4 1st Qu.:3.504 1st Qu.:1.6464 1st Qu.:1.767
## Median :299.9 Median :3.620 Median :1.7411 Median :1.853
## Mean :299.9 Mean :3.620 Mean :1.6833 Mean :1.937
## 3rd Qu.:300.3 3rd Qu.:3.701 3rd Qu.:1.7992 3rd Qu.:2.012
## Max. :302.4 Max. :4.194 Max. :1.9326 Max. :3.021
## NA's :536 NA's :536 NA's :536 NA's :536
## SMI temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
## Min. : 7.317 Min. :23.30 Min. :0.5966 Min. :13.75
## 1st Qu.: 9.748 1st Qu.:24.05 1st Qu.:0.7828 1st Qu.:29.21
## Median :10.624 Median :29.74 Median :1.0461 Median :45.24
## Mean :10.607 Mean :29.61 Mean :1.1502 Mean :52.95
## 3rd Qu.:11.348 3rd Qu.:35.05 3rd Qu.:1.5191 3rd Qu.:82.84
## Max. :14.263 Max. :36.00 Max. :1.8447 Max. :93.15
##
## humidity_SD_tmttrial e_a_mean_tmttrial e_a_SD_tmttrial VPD_mean_tmttrial
## Min. : 4.370 Min. :0.4258 Min. :0.1891 Min. :0.2053
## 1st Qu.: 6.234 1st Qu.:1.3896 1st Qu.:0.2651 1st Qu.:0.8018
## Median : 7.382 Median :2.2917 Median :0.3166 Median :2.0450
## Mean : 8.758 Mean :2.3359 Mean :0.3934 Mean :2.0051
## 3rd Qu.:12.297 3rd Qu.:2.6749 3rd Qu.:0.5019 3rd Qu.:3.1651
## Max. :19.846 Max. :4.8058 Max. :0.7691 Max. :4.0685
##
## VPD_SD_tmttrial temp_C humidity_percent e_a_kPa
## Min. :0.1690 Min. :23.80 Min. :17.90 Min. :0.500
## 1st Qu.:0.2558 1st Qu.:23.80 1st Qu.:34.10 1st Qu.:2.000
## Median :0.3483 Median :29.65 Median :56.20 Median :2.150
## Mean :0.4036 Mean :29.65 Mean :52.85 Mean :2.322
## 3rd Qu.:0.4735 3rd Qu.:35.50 3rd Qu.:78.30 3rd Qu.:2.300
## Max. :0.8863 Max. :35.50 Max. :80.90 Max. :4.500
##
## VPD_kPa
## Min. :0.600
## 1st Qu.:0.600
## Median :1.800
## Mean :2.003
## 3rd Qu.:3.800
## Max. :3.800
## dat_no_rehab %>%
group_by(type) %>%
summarise(mean(msmt_temp_C, na.rm = T),
sd(msmt_temp_C, na.rm = T),
mean(msmt_RH_percent, na.rm = T),
sd(msmt_RH_percent, na.rm = T),
mean(msmt_VPD_kPa, na.rm = T),
sd(msmt_VPD_kPa, na.rm = T))## # A tibble: 1 × 7
## type `mean(msmt_temp_C, na.rm = T)` sd(msmt…¹ mean(…² sd(ms…³ mean(…⁴ sd(ms…⁵
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 exp 26.7 0.799 46.7 6.76 1.94 0.337
## # … with abbreviated variable names ¹`sd(msmt_temp_C, na.rm = T)`,
## # ²`mean(msmt_RH_percent, na.rm = T)`, ³`sd(msmt_RH_percent, na.rm = T)`,
## # ⁴`mean(msmt_VPD_kPa, na.rm = T)`, ⁵`sd(msmt_VPD_kPa, na.rm = T)`Means by Day
Calculate mean values per day per tmt group.
means <- dat_no_rehab %>%
group_by(day_n, tmt) %>%
summarise(n_lizards = n(),
mean_CEWL = mean(CEWL_g_m2h_mean, na.rm = T),
sd_CEWL = sd(CEWL_g_m2h_mean, na.rm = T),
mean_osml = mean(osmolality_mmol_kg_mean, na.rm = T),
sd_osml = sd(osmolality_mmol_kg_mean, na.rm = T),
mean_hct = mean(hematocrit_percent, na.rm = T),
sd_hct = sd(hematocrit_percent, na.rm = T),
mean_mass = mean(mass_g, na.rm = T),
sd_mass = sd(mass_g, na.rm = T),
mean_SMI = mean(SMI, na.rm = T),
sd_SMI = sd(SMI, na.rm = T)) %>%
mutate(se_CEWL = (sd_CEWL/sqrt(n_lizards)),
se_osml = (sd_osml/sqrt(n_lizards)),
se_hct = (sd_hct/sqrt(n_lizards)),
se_mass = (sd_mass/sqrt(n_lizards)),
se_SMI = (sd_SMI/sqrt(n_lizards)))## `summarise()` has grouped output by 'day_n'. You can override using the
## `.groups` argument.# get rid of non-defined points
# these are from when not all variables were measured for a given date
means$mean_CEWL[is.nan(means$mean_CEWL)] <- NA
means$mean_osml[is.nan(means$mean_osml)] <- NA
means$mean_hct[is.nan(means$mean_hct)] <- NA
means$mean_mass[is.nan(means$mean_mass)] <- NA
means$mean_SMI[is.nan(means$mean_SMI)] <- NA
means## # A tibble: 24 × 18
## # Groups: day_n [6]
## day_n tmt n_liz…¹ mean_…² sd_CEWL mean_…³ sd_osml mean_…⁴ sd_hct mean_…⁵
## <dbl> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 Cool Hu… 34 20.9 4.78 351. 20.3 39.6 5.30 11.6
## 2 0 Hot Hum… 33 21.4 4.85 347. 18.7 37.9 5.46 11.6
## 3 0 Cool Dry 33 20.0 6.08 346. 20.6 39.3 5.96 11.8
## 4 0 Hot Dry 34 20.9 5.93 347. 17.9 38.9 5.05 12.0
## 5 4 Cool Hu… 34 NA NA 356. 24.5 34.5 5.34 11.1
## 6 4 Hot Hum… 33 NA NA 359. 22.5 32.0 5.38 10.6
## 7 4 Cool Dry 33 NA NA 355. 27.0 35.0 7.02 11.1
## 8 4 Hot Dry 34 NA NA 361. 25.8 33.1 5.00 10.6
## 9 5 Cool Hu… 34 NA NA NA NA NA NA 10.8
## 10 5 Hot Hum… 33 NA NA NA NA NA NA 10.2
## # … with 14 more rows, 8 more variables: sd_mass <dbl>, mean_SMI <dbl>,
## # sd_SMI <dbl>, se_CEWL <dbl>, se_osml <dbl>, se_hct <dbl>, se_mass <dbl>,
## # se_SMI <dbl>, and abbreviated variable names ¹n_lizards, ²mean_CEWL,
## # ³mean_osml, ⁴mean_hct, ⁵mean_massEnd Values Only
Select all raw measurements for only day=8 values.
delta CEWL
Get a df that only has complete observations that include CEWL values (only obs from before and after the experiment). Then, calculate the change (delta) in CEWL from before to after the experiment. Because we only measured CEWL at those two time points, it makes more sense to assess the amount of change in CEWL for each lizard, rather than measuring the change over time.
start_CEWL <- dat_no_rehab %>%
dplyr::filter(day_n == 0) %>%
dplyr::select(individual_ID, start_CEWL = CEWL_g_m2h_mean)
dat_no_rehab_deltaCEWL <- dat_no_rehab %>% # initiate new df
dplyr::filter(complete.cases(CEWL_g_m2h_mean)) %>% # only use obs incl CEWL
dplyr::filter(day_n == 8) %>% # get only obs for post-exp
left_join(start_CEWL, by = 'individual_ID') %>% # add start CEWL to both obs for each lizard
mutate(delta_CEWL = CEWL_g_m2h_mean - start_CEWL) # calculate deltaCEWL after-before experimentExperiment Models
We predicted that there would be effects of day, humidity treatment, temperature treatment, and treatment VPD. However, we can’t use the standard backwards model selection because the three treatment variables are collinear (VIF much higher than acceptable) and it leads to issues with changing the sign of the estimates when all three are included together. So, we will run singular models with each treatment variable alone:
response ~ dayhumidity response ~ daytemperature response ~ day*VPD
I added models with water vapor pressure based on comments received during review, but deemed it uninformative to add to the results presented in the paper
Then, we will use AIC, RMSE, and R-sq to assess which treatment effect is most important to that response variable.
Body Mass
Building
Build each treatment effect model.
mass_mod_WVP <- lmerTest::lmer(data = dat_no_rehab,
mass_g ~ day_n * e_a_kPa +
(1|trial_number/individual_ID))## boundary (singular) fit: see help('isSingular')mass_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
mass_g ~ day_n * VPD_kPa +
(1|trial_number/individual_ID))## boundary (singular) fit: see help('isSingular')## Warning: Model failed to converge with 1 negative eigenvalue: -1.7e+01mass_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
mass_g ~ day_n * humidity_tmt +
(1|trial_number/individual_ID))## boundary (singular) fit: see help('isSingular')mass_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
mass_g ~ day_n * temp_tmt +
(1|trial_number/individual_ID))## boundary (singular) fit: see help('isSingular')Assumptions
Check linear regression assumptions/conditions.
##
## Shapiro-Wilk normality test
##
## data: residuals(mass_mod_VPD)
## W = 0.93331, p-value < 2.2e-16
##
## Shapiro-Wilk normality test
##
## data: residuals(mass_mod_hum)
## W = 0.93039, p-value < 2.2e-16
##
## Shapiro-Wilk normality test
##
## data: residuals(mass_mod_temp)
## W = 0.90802, p-value < 2.2e-16Normality is violated, but linearity, equal error variance, and independence are all good.
Comparison
Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.
We…
- calculate RMSE manually
- use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
- use the aictab function in the AICmodavg package to get AIC and deltaAIC values
- get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
mass_RMSE_Rsq <- data.frame(model =
c('Day * WVP',
'Day * VPD',
'Day * Humidity',
'Day * Temp'
),
RMSE = c(sqrt(mean((residuals(mass_mod_WVP))^2)),
sqrt(mean((residuals(mass_mod_VPD))^2)),
sqrt(mean((residuals(mass_mod_hum))^2)),
sqrt(mean((residuals(mass_mod_temp))^2))),
# marginal Rsq for the amount of variance
# explained by fixed effects only
Rsq = c(MuMIn::r.squaredGLMM(mass_mod_WVP)[,"R2m"],
MuMIn::r.squaredGLMM(mass_mod_VPD)[,"R2m"],
MuMIn::r.squaredGLMM(mass_mod_hum)[,"R2m"],
MuMIn::r.squaredGLMM(mass_mod_temp)[,"R2m"]))## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.# calculate AIC
mass_models <- list(mass_mod_WVP, mass_mod_VPD, mass_mod_hum, mass_mod_temp)
EXP_mod_names <- data.frame(model =
c('Day * WVP',
'Day * VPD',
'Day * Humidity',
'Day * Temp'
))
mass_AICc <- data.frame(aictab(cand.set = mass_models,
modnames = EXP_mod_names))## Warning in aictab.AIClmerModLmerTest(cand.set = mass_models, modnames = EXP_mod_names):
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects# compare across models
mass_across <- mass_RMSE_Rsq %>%
left_join(mass_AICc, by = 'model') %>%
mutate(response = "Body Mass (g)") %>%
arrange(Delta_AICc)
# calculate F & p-values
mass_VPD_p <- data.frame(anova(mass_mod_VPD,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * VPD',
predictor = rownames(.))
mass_hum_p <- data.frame(anova(mass_mod_hum,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Humidity',
predictor = rownames(.))
mass_temp_p <- data.frame(anova(mass_mod_temp,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Temp',
predictor = rownames(.))
# save within model values
mass_within <- mass_VPD_p %>%
rbind(mass_hum_p) %>%
rbind(mass_temp_p) %>%
mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
response = "Body Mass (g)")Body Condition
not in publication
Building
Build each treatment effect model.
SMI_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
SMI ~ day_n * VPD_kPa +
(1|trial_number/individual_ID))
SMI_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
SMI ~ day_n * humidity_tmt +
(1|trial_number/individual_ID))
SMI_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
SMI ~ day_n * temp_tmt +
(1|trial_number/individual_ID))Assumptions
Check linear regression assumptions/conditions.
##
## Shapiro-Wilk normality test
##
## data: residuals(SMI_mod_VPD)
## W = 0.96044, p-value = 6.428e-14
##
## Shapiro-Wilk normality test
##
## data: residuals(SMI_mod_hum)
## W = 0.95569, p-value = 7.603e-15
##
## Shapiro-Wilk normality test
##
## data: residuals(SMI_mod_temp)
## W = 0.93877, p-value < 2.2e-16Normality is violated, but linearity, equal error variance, and independence are all good.
Comparison
Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.
We…
- calculate RMSE manually
- use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
- use the aictab function in the AICmodavg package to get AIC and deltaAIC values
- get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
SMI_RMSE_Rsq <- data.frame(model =
c('Day * VPD',
'Day * Humidity',
'Day * Temp'
),
RMSE = c(sqrt(mean((residuals(SMI_mod_VPD))^2)),
sqrt(mean((residuals(SMI_mod_hum))^2)),
sqrt(mean((residuals(SMI_mod_temp))^2))),
# marginal Rsq for the amount of variance
# explained by fixed effects only
Rsq = c(MuMIn::r.squaredGLMM(SMI_mod_VPD)[,"R2m"],
MuMIn::r.squaredGLMM(SMI_mod_hum)[,"R2m"],
MuMIn::r.squaredGLMM(SMI_mod_temp)[,"R2m"]))
# calculate AIC
SMI_models <- list(SMI_mod_VPD, SMI_mod_hum, SMI_mod_temp)
EXP_mod_names <- data.frame(model =
c('Day * VPD',
'Day * Humidity',
'Day * Temp'
))
SMI_AICc <- data.frame(aictab(cand.set = SMI_models,
modnames = EXP_mod_names))## Warning in aictab.AIClmerModLmerTest(cand.set = SMI_models, modnames = EXP_mod_names):
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects# compare across models
SMI_across <- SMI_RMSE_Rsq %>%
left_join(SMI_AICc, by = 'model') %>%
mutate(response = "Body Condition (g')") %>%
arrange(Delta_AICc)
# calculate F & p-values
SMI_VPD_p <- data.frame(anova(SMI_mod_VPD,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * VPD',
predictor = rownames(.))
SMI_hum_p <- data.frame(anova(SMI_mod_hum,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Humidity',
predictor = rownames(.))
SMI_temp_p <- data.frame(anova(SMI_mod_temp,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Temp',
predictor = rownames(.))
# save within model values
SMI_within <- SMI_VPD_p %>%
rbind(SMI_hum_p) %>%
rbind(SMI_temp_p) %>%
mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
response = "Body Condition (g')")Hematocrit
Building
Build each treatment effect model.
hct_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
hematocrit_percent ~ day_n * VPD_kPa +
(1|trial_number/individual_ID))
hct_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
hematocrit_percent ~ day_n * humidity_tmt +
(1|trial_number/individual_ID))
hct_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
hematocrit_percent ~ day_n * temp_tmt +
(1|trial_number/individual_ID))Assumptions
Check linear regression assumptions/conditions.
##
## Shapiro-Wilk normality test
##
## data: residuals(hct_mod_VPD)
## W = 0.99652, p-value = 0.5454
##
## Shapiro-Wilk normality test
##
## data: residuals(hct_mod_hum)
## W = 0.99619, p-value = 0.4584
##
## Shapiro-Wilk normality test
##
## data: residuals(hct_mod_temp)
## W = 0.99478, p-value = 0.1975All assumptions, normality, linearity, equal error variance, and independence are all good.
Comparison
Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.
# calculate RMSE & R^2
hct_RMSE_Rsq <- data.frame(model =
c('Day * VPD',
'Day * Humidity',
'Day * Temp'
),
RMSE = c(sqrt(mean((residuals(hct_mod_VPD))^2)),
sqrt(mean((residuals(hct_mod_hum))^2)),
sqrt(mean((residuals(hct_mod_temp))^2))),
# marginal Rsq for the amount of variance
# explained by fixed effects only
Rsq = c(MuMIn::r.squaredGLMM(hct_mod_VPD)[,"R2m"],
MuMIn::r.squaredGLMM(hct_mod_hum)[,"R2m"],
MuMIn::r.squaredGLMM(hct_mod_temp)[,"R2m"]))
# calculate AIC
hct_models <- list(hct_mod_VPD, hct_mod_hum, hct_mod_temp)
hct_AICc <- data.frame(aictab(cand.set = hct_models,
modnames = EXP_mod_names))## Warning in aictab.AIClmerModLmerTest(cand.set = hct_models, modnames = EXP_mod_names):
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects# compare across models
hct_across <- hct_RMSE_Rsq %>%
left_join(hct_AICc, by = 'model') %>%
mutate(response = "Hematocrit (%)") %>%
arrange(Delta_AICc)
# calculate F & p-values
hct_VPD_p <- data.frame(anova(hct_mod_VPD,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * VPD',
predictor = rownames(.))
hct_hum_p <- data.frame(anova(hct_mod_hum,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Humidity',
predictor = rownames(.))
hct_temp_p <- data.frame(anova(hct_mod_temp,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Temp',
predictor = rownames(.))
# save within model values
hct_within <- hct_VPD_p %>%
rbind(hct_hum_p) %>%
rbind(hct_temp_p) %>%
mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
response = "Hematocrit (%)")Osmolality
Building
Build each treatment effect model.
osml_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
osmolality_mmol_kg_mean ~ day_n * VPD_kPa +
(1|trial_number/individual_ID))
osml_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
osmolality_mmol_kg_mean ~ day_n * humidity_tmt +
(1|trial_number/individual_ID))
osml_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
osmolality_mmol_kg_mean ~ day_n * temp_tmt +
(1|trial_number/individual_ID))Assumptions
Check linear regression assumptions/conditions.
##
## Shapiro-Wilk normality test
##
## data: residuals(osml_mod_VPD)
## W = 0.97687, p-value = 6.755e-06
##
## Shapiro-Wilk normality test
##
## data: residuals(osml_mod_hum)
## W = 0.97746, p-value = 8.914e-06
##
## Shapiro-Wilk normality test
##
## data: residuals(osml_mod_temp)
## W = 0.97346, p-value = 1.44e-06Normality is violated, but linearity, equal error variance, and independence are all okay.
Comparison
Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.
# calculate RMSE & R^2
osml_RMSE_Rsq <- data.frame(model =
c('Day * VPD',
'Day * Humidity',
'Day * Temp'
),
RMSE = c(sqrt(mean((residuals(osml_mod_VPD))^2)),
sqrt(mean((residuals(osml_mod_hum))^2)),
sqrt(mean((residuals(osml_mod_temp))^2))),
# marginal Rsq for the amount of variance
# explained by fixed effects only
Rsq = c(MuMIn::r.squaredGLMM(osml_mod_VPD)[,"R2m"],
MuMIn::r.squaredGLMM(osml_mod_hum)[,"R2m"],
MuMIn::r.squaredGLMM(osml_mod_temp)[,"R2m"]))
# calculate AIC
osml_models <- list(osml_mod_VPD, osml_mod_hum, osml_mod_temp)
osml_AICc <- data.frame(aictab(cand.set = osml_models,
modnames = EXP_mod_names))## Warning in aictab.AIClmerModLmerTest(cand.set = osml_models, modnames = EXP_mod_names):
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects# compare across models
osml_across <- osml_RMSE_Rsq %>%
left_join(osml_AICc, by = 'model') %>%
mutate(response = "Plasma Osmolality (mmol/kg)") %>%
arrange(Delta_AICc)
# calculate F & p-values
osml_VPD_p <- data.frame(anova(osml_mod_VPD,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * VPD',
predictor = rownames(.))
osml_hum_p <- data.frame(anova(osml_mod_hum,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Humidity',
predictor = rownames(.))
osml_temp_p <- data.frame(anova(osml_mod_temp,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Day * Temp',
predictor = rownames(.))
# save within model values
osml_within <- osml_VPD_p %>%
rbind(osml_hum_p) %>%
rbind(osml_temp_p) %>%
mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
response = "Plasma Osmolality (mmol/kg)")Body Temp
I need to double check whether CEWL has a relationship with body temperature at the point of measurement.
## Warning: Removed 537 rows containing missing values (`geom_point()`).
Test an lm of raw CEWL ~ body temp for day 8 measurements only.
body_temp_test_dat <- end_vals %>%
dplyr::filter(complete.cases(CEWL_g_m2h_mean))
CEWL_body_temp_mod <- lmerTest::lmer(data = body_temp_test_dat,
CEWL_g_m2h_mean ~ cloacal_temp_C * tmt +
(1|trial_number))
summary(CEWL_body_temp_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ cloacal_temp_C * tmt + (1 | trial_number)
## Data: body_temp_test_dat
##
## REML criterion at convergence: 839.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.02426 -0.58867 -0.09161 0.46172 3.09926
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 4.977 2.231
## Residual 36.926 6.077
## Number of obs: 133, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 42.83759 26.65061 123.21107 1.607 0.111
## cloacal_temp_C -0.49420 1.05450 122.93003 -0.469 0.640
## tmtHot Humid -16.46275 34.17703 122.09739 -0.482 0.631
## tmtCool Dry -23.23095 47.80497 123.10602 -0.486 0.628
## tmtHot Dry -3.95989 36.73317 123.70420 -0.108 0.914
## cloacal_temp_C:tmtHot Humid 0.90966 1.35821 122.13178 0.670 0.504
## cloacal_temp_C:tmtCool Dry 0.65396 1.86629 123.06213 0.350 0.727
## cloacal_temp_C:tmtHot Dry -0.06271 1.44351 123.65911 -0.043 0.965
##
## Correlation of Fixed Effects:
## (Intr) clc__C tmtHtH tmtClD tmtHtD c__C:H c__C:CD
## clocl_tmp_C -0.999
## tmtHotHumid -0.745 0.744
## tmtCool Dry -0.528 0.528 0.434
## tmtHot Dry -0.694 0.694 0.564 0.422
## clcl_t_C:HH 0.741 -0.741 -0.999 -0.432 -0.562
## clcl_t_C:CD 0.536 -0.536 -0.439 -0.999 -0.427 0.438
## clcl_t_C:HD 0.699 -0.700 -0.567 -0.425 -0.999 0.566 0.430## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## cloacal_temp_C 1.1914 1.1914 1 115.35 0.0323 0.8578
## tmt 14.9005 4.9668 3 122.58 0.1345 0.9393
## cloacal_temp_C:tmt 27.5730 9.1910 3 122.56 0.2489 0.8620ggplot(body_temp_test_dat) +
geom_point(aes(x = msmt_temp_C,
y = cloacal_temp_C)) +
xlim(23, 28) + ylim(23,28)
Do the same exact stats for ambient temp and VPD at the time of measurement:
CEWL_msmt_temp_mod <- lmerTest::lmer(data = body_temp_test_dat,
CEWL_g_m2h_mean ~ msmt_temp_C * tmt +
(1|trial_number))
CEWL_msmt_WVP_mod <- lmerTest::lmer(data = body_temp_test_dat,
CEWL_g_m2h_mean ~ e_a_kPa_m * tmt +
(1|trial_number))
CEWL_msmt_VPD_mod <- lmerTest::lmer(data = body_temp_test_dat,
CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt +
(1|trial_number))
summary(CEWL_msmt_temp_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_temp_C * tmt + (1 | trial_number)
## Data: body_temp_test_dat
##
## REML criterion at convergence: 824.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8115 -0.5673 -0.0847 0.4720 3.5376
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 7.298 2.701
## Residual 33.499 5.788
## Number of obs: 133, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -39.6305 48.2297 112.4351 -0.822 0.413
## msmt_temp_C 2.6733 1.8380 112.5242 1.454 0.149
## tmtHot Humid -73.4611 65.2521 121.4757 -1.126 0.262
## tmtCool Dry -18.6586 60.6146 121.3298 -0.308 0.759
## tmtHot Dry -61.0925 64.0861 121.2408 -0.953 0.342
## msmt_temp_C:tmtHot Humid 3.0314 2.4837 121.4794 1.221 0.225
## msmt_temp_C:tmtCool Dry 0.4532 2.3073 121.3153 0.196 0.845
## msmt_temp_C:tmtHot Dry 2.1152 2.4428 121.2415 0.866 0.388
##
## Correlation of Fixed Effects:
## (Intr) msm__C tmtHtH tmtClD tmtHtD m__C:H m__C:CD
## msmt_temp_C -0.999
## tmtHotHumid -0.593 0.592
## tmtCool Dry -0.614 0.613 0.467
## tmtHot Dry -0.580 0.580 0.445 0.478
## msmt_t_C:HH 0.594 -0.594 -1.000 -0.467 -0.445
## msmt_t_C:CD 0.614 -0.614 -0.467 -1.000 -0.478 0.467
## msmt_t_C:HD 0.579 -0.579 -0.444 -0.478 -1.000 0.445 0.478## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt + (1 | trial_number)
## Data: body_temp_test_dat
##
## REML criterion at convergence: 817.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8187 -0.5552 -0.0438 0.3596 3.2064
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 11.75 3.428
## Residual 34.18 5.846
## Number of obs: 133, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.583 14.417 63.780 1.497 0.1393
## msmt_VPD_kPa 4.952 8.054 67.896 0.615 0.5407
## tmtHot Humid -24.082 16.730 120.618 -1.439 0.1526
## tmtCool Dry -9.284 16.225 120.721 -0.572 0.5682
## tmtHot Dry -25.308 16.211 120.503 -1.561 0.1211
## msmt_VPD_kPa:tmtHot Humid 17.094 9.373 120.640 1.824 0.0707 .
## msmt_VPD_kPa:tmtCool Dry 1.504 9.025 120.677 0.167 0.8679
## msmt_VPD_kPa:tmtHot Dry 11.097 9.093 120.529 1.220 0.2247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) m_VPD_ tmtHtH tmtClD tmtHtD m_VPDH m_VPD_P:CD
## msmt_VPD_kP -0.992
## tmtHotHumid -0.538 0.538
## tmtCool Dry -0.524 0.523 0.472
## tmtHot Dry -0.535 0.535 0.476 0.490
## ms_VPD_P:HH 0.538 -0.542 -0.996 -0.470 -0.474
## ms_VPD_P:CD 0.527 -0.531 -0.474 -0.996 -0.492 0.476
## ms_VPD_P:HD 0.529 -0.533 -0.473 -0.488 -0.996 0.475 0.493## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## msmt_temp_C 309.546 309.546 1 59.634 9.2405 0.003513 **
## tmt 57.484 19.161 3 121.536 0.5720 0.634497
## msmt_temp_C:tmt 65.354 21.785 3 121.538 0.6503 0.584244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## e_a_kPa_m 141.13 141.126 1 8.973 4.0770 0.07432 .
## tmt 174.69 58.230 3 121.755 1.6822 0.17438
## e_a_kPa_m:tmt 129.98 43.325 3 121.753 1.2516 0.29414
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## msmt_VPD_kPa 115.03 115.026 1 32.112 3.3653 0.07586 .
## tmt 112.93 37.643 3 121.388 1.1013 0.35140
## msmt_VPD_kPa:tmt 152.90 50.966 3 121.397 1.4911 0.22039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL
First, get a separate df of end-experiment delta CEWL for each tmt group:
dat_no_rehab_deltaCEWL_post <- dat_no_rehab_deltaCEWL %>%
dplyr::filter(day_n == 8)
HH8 <- dat_no_rehab_deltaCEWL_post %>%
dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD8 <- dat_no_rehab_deltaCEWL_post %>%
dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH8 <- dat_no_rehab_deltaCEWL_post %>%
dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD8 <- dat_no_rehab_deltaCEWL_post %>%
dplyr::filter(substr(tmt, 1, 6) == "Cool D")t-tests
I use t-tests to assess whether the change in CEWL from the experiment is significantly different from zero for each treatment group.
##
## One Sample t-test
##
## data: HH8$delta_CEWL
## t = 8.7341, df = 32, p-value = 5.573e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 11.76821 18.92672
## sample estimates:
## mean of x
## 15.34746##
## One Sample t-test
##
## data: HD8$delta_CEWL
## t = 2.513, df = 33, p-value = 0.01703
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.7042497 6.6929758
## sample estimates:
## mean of x
## 3.698613##
## One Sample t-test
##
## data: CH8$delta_CEWL
## t = 8.7488, df = 32, p-value = 5.363e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 7.296566 11.725292
## sample estimates:
## mean of x
## 9.510929##
## One Sample t-test
##
## data: CD8$delta_CEWL
## t = 2.75, df = 32, p-value = 0.009721
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.9302207 6.2446177
## sample estimates:
## mean of x
## 3.587419Building
Build each treatment effect model. These are equivalent to the models built for body condition, hematocrit, and plasma osmolality, but the effect of time is not included because we assess CEWL as delta CEWL.
CEWL_mod_VPD <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ VPD_kPa +
(1|trial_number))
CEWL_mod_WVP <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ e_a_kPa +
(1|trial_number))
CEWL_mod_hum <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ humidity_tmt +
(1|trial_number))
CEWL_mod_temp <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ temp_tmt +
(1|trial_number))
CEWL_mod_double <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ temp_tmt*humidity_tmt +
(1|trial_number))
# summary(CEWL_mod_VPD)
# summary(CEWL_mod_hum)
# summary(CEWL_mod_temp)
# summary(CEWL_mod_double)Assumptions
Check linear regression assumptions/conditions.
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_mod_VPD)
## W = 0.97248, p-value = 0.008418
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_mod_hum)
## W = 0.98319, p-value = 0.1002
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_mod_temp)
## W = 0.98262, p-value = 0.08761All assumptions are fine.
Comparison
Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.
# calculate RMSE & R^2
CEWL_RMSE_Rsq <- data.frame(model =
c('WVP',
'VPD',
'Humidity',
'Temp'
),
RMSE = c(sqrt(mean((residuals(CEWL_mod_WVP))^2)),
sqrt(mean((residuals(CEWL_mod_VPD))^2)),
sqrt(mean((residuals(CEWL_mod_hum))^2)),
sqrt(mean((residuals(CEWL_mod_temp))^2))),
# marginal Rsq for the amount of variance
# explained by fixed effects only
Rsq = c(MuMIn::r.squaredGLMM(CEWL_mod_WVP)[,"R2m"],
MuMIn::r.squaredGLMM(CEWL_mod_VPD)[,"R2m"],
MuMIn::r.squaredGLMM(CEWL_mod_hum)[,"R2m"],
MuMIn::r.squaredGLMM(CEWL_mod_temp)[,"R2m"]))
# calculate AIC
CEWL_models <- list(CEWL_mod_WVP, CEWL_mod_VPD, CEWL_mod_hum, CEWL_mod_temp)
CEWL_mod_names <- data.frame(model =
c('WVP',
'VPD',
'Humidity',
'Temp'
))
CEWL_AICc <- data.frame(aictab(cand.set = CEWL_models,
modnames = CEWL_mod_names))## Warning in aictab.AIClmerModLmerTest(cand.set = CEWL_models, modnames = CEWL_mod_names):
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects# compare across models
CEWL_across <- CEWL_RMSE_Rsq %>%
left_join(CEWL_AICc, by = 'model') %>%
mutate(response = "deltaCEWL") %>%
arrange(Delta_AICc)
# calculate F & p-values
CEWL_VPD_p <- data.frame(anova(CEWL_mod_VPD,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'VPD',
predictor = rownames(.))
CEWL_hum_p <- data.frame(anova(CEWL_mod_hum,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Humidity',
predictor = rownames(.))
CEWL_temp_p <- data.frame(anova(CEWL_mod_temp,
type = "3",
ddf = "Kenward-Roger")) %>%
mutate(model = 'Temp',
predictor = rownames(.))
# save within model values
CEWL_within <- CEWL_VPD_p %>%
rbind(CEWL_hum_p) %>%
rbind(CEWL_temp_p) %>%
mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
response = "deltaCEWL")But, water vapor pressure is pretty confounded with humidity tmt.
Group Export
Put all the model statistics into one df/csv - one for among-model comparisons, and one for within-model stats.
# comparisons of all models for each variable
experiment_model_compare <- CEWL_across %>%
rbind(osml_across) %>%
rbind(hct_across) %>%
rbind(mass_across) %>%
rbind(SMI_across) %>%
dplyr::select(response, model,
RMSE, Rsq, AICc, Delta_AICc) %>%
mutate(RMSE = round(RMSE, 2),
Rsq = round(Rsq, 2),
AICc = round(AICc, 2),
Delta_AICc = round(Delta_AICc, 2))
write.csv(experiment_model_compare,
"./results_statistics/exp_model_comparisons.csv")
# all of the actual model results
experiment_model_values <- CEWL_within %>%
rbind(osml_within) %>%
rbind(hct_within) %>%
rbind(mass_within) %>%
rbind(SMI_within) %>%
dplyr::select(response, model, predictor,
seq_sum_of_squares = Sum.Sq,
df,
F_statistic = F.value,
p_value = Pr..F.) %>%
mutate(seq_sum_of_squares = round(seq_sum_of_squares, 0),
F_statistic = round(F_statistic, 2))
write.csv(experiment_model_values, "./results_statistics/exp_model_values.csv")Effect Estimates
End Value CIs
Now, we can use the emmeans function from the emmeans package to estimate marginal means and confidence intervals for the values among treatment groups at the end of the experiment. But, to get for each treatment group, we need to run a new model with day 8 data only and tmt as a single, 4-category variable.
## boundary (singular) fit: see help('isSingular')mass_emmeans <- data.frame(emmeans(mass_mod_end, "tmt")) %>%
mutate(response = "Body Mass (g)")
mass_pairwise <- data.frame(pairs(emmeans(mass_mod_end, "tmt"))) %>%
mutate(response = "Body Mass (g)")
# Body Condition
SMI_mod_end <- lmerTest::lmer(data = end_vals,
SMI ~ tmt +
(1|trial_number))
SMI_emmeans <- data.frame(emmeans(SMI_mod_end, "tmt")) %>%
mutate(response = "Body Condition (g')")
SMI_pairwise <- data.frame(pairs(emmeans(SMI_mod_end, "tmt"))) %>%
mutate(response = "Body Condition (g')")
# Hematocrit
hct_mod_end <- lmerTest::lmer(data = end_vals,
hematocrit_percent ~ tmt +
(1|trial_number))
hct_emmeans <- data.frame(emmeans(hct_mod_end, "tmt")) %>%
mutate(response = "Hematocrit (%)")
hct_pairwise <- data.frame(pairs(emmeans(hct_mod_end, "tmt"))) %>%
mutate(response = "Hematocrit (%)")
# Plasma Osmolality
osml_mod_end <- lmerTest::lmer(data = end_vals,
osmolality_mmol_kg_mean ~ tmt +
(1|trial_number))
osml_emmeans <- data.frame(emmeans(osml_mod_end, "tmt")) %>%
mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairwise <- data.frame(pairs(emmeans(osml_mod_end, "tmt"))) %>%
mutate(response = "Plasma Osmolality (mmol/kg)")
# CEWL
CEWL_mod_end <- lmerTest::lmer(data = end_vals,
CEWL_g_m2h_mean ~ tmt +
(1|trial_number))
CEWL_emmeans <- data.frame(emmeans(CEWL_mod_end, "tmt")) %>%
mutate(response = "CEWL (g/m2h)")
CEWL_pairwise <- data.frame(pairs(emmeans(CEWL_mod_end, "tmt"))) %>%
mutate(response = "CEWL (g/m2h)")Rate of Change
# Body MASS
mass_mod_day <- lmerTest::lmer(data = dat_no_rehab,
mass_g ~ day_n * tmt +
(1|trial_number/individual_ID))## boundary (singular) fit: see help('isSingular')mass_emtrends <- data.frame(emtrends(mass_mod_day, "tmt", var = "day_n")) %>%
mutate(response = "Body Mass (g)")
mass_pairtrend <- data.frame(pairs(emtrends(mass_mod_day, "tmt", var = "day_n"))) %>%
mutate(response = "Body Mass (g)")
# Body Condition
SMI_mod_day <- lmerTest::lmer(data = dat_no_rehab,
SMI ~ day_n * tmt +
(1|trial_number/individual_ID))
SMI_emtrends <- data.frame(emtrends(SMI_mod_day, "tmt", var = "day_n")) %>%
mutate(response = "Body Condition (g')")
SMI_pairtrend <- data.frame(pairs(emtrends(SMI_mod_day, "tmt", var = "day_n"))) %>%
mutate(response = "Body Condition (g')")
# Hematocrit
hct_mod_day <- lmerTest::lmer(data = dat_no_rehab,
hematocrit_percent ~ day_n * tmt +
(1|trial_number/individual_ID))
hct_emtrends <- data.frame(emtrends(hct_mod_day, "tmt", var = "day_n")) %>%
mutate(response = "Hematocrit (%)")
hct_pairtrend <- data.frame(pairs(emtrends(hct_mod_day, "tmt", var = "day_n"))) %>%
mutate(response = "Hematocrit (%)")
# Plasma Osmolality
osml_mod_day <- lmerTest::lmer(data = dat_no_rehab,
osmolality_mmol_kg_mean ~ day_n * tmt +
(1|trial_number/individual_ID))
osml_emtrends <- data.frame(emtrends(osml_mod_day, "tmt", var = "day_n")) %>%
mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairtrend <- data.frame(pairs(emtrends(osml_mod_day, "tmt", var = "day_n"))) %>%
mutate(response = "Plasma Osmolality (mmol/kg)")
# CEWL
CEWL_mod_day <- lmerTest::lmer(data = dat_no_rehab,
CEWL_g_m2h_mean ~ day_n * tmt +
(1|trial_number))
CEWL_emtrends <- data.frame(emtrends(CEWL_mod_day, "tmt", var = "day_n")) %>%
mutate(response = "CEWL (g/m2h)")
CEWL_pairtrend <- data.frame(pairs(emtrends(CEWL_mod_day, "tmt", var = "day_n"))) %>%
mutate(response = "CEWL (g/m2h)")
# put together
all_emtrends <- CEWL_emtrends %>%
rbind(osml_emtrends) %>%
rbind(hct_emtrends) %>%
rbind(mass_emtrends) %>%
rbind(SMI_emtrends) %>%
mutate(confint95 = paste(round(lower.CL, 2), round(upper.CL, 2), sep = ", ")) %>%
dplyr::select(response, tmt,
per_day_change = day_n.trend,
confint95,
SE, df)
write.csv(all_emtrends, "./results_statistics/exp_emtrends_per_day_change.csv")Given the daily trends for plasma osmolality, in total for acclimation change in osml had 1.04 (0.13x8) to 13 (1.587x8) mmol/kg of change.
CEWL ~ osmolality
Let’s compare the relationship of CEWL ~ osmolality post-experiment to what we found pre-experiment in the “analysis_capture” file. The treatment-specific sub-dataframes were created in the “Body Temp” model subsection above.
These are NOT the stats presents for Figure 6
CEWL_osml <- lmerTest::lmer(data = end_vals,
CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
(1|trial_number))
summary(CEWL_osml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
## Data: end_vals
##
## REML criterion at convergence: 855.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1565 -0.6405 -0.2787 0.6163 3.5785
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 5.66 2.379
## Residual 58.25 7.632
## Number of obs: 123, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 18.53375 10.95312 66.51336 1.692 0.0953 .
## osmolality_mmol_kg_mean 0.02904 0.03081 72.84361 0.942 0.3491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.993## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 45.195 45.195 1 69.109 0.7759 0.3814HH_CEWL_osml <- lmerTest::lmer(data = HH8,
CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
(1|trial_number))## boundary (singular) fit: see help('isSingular')## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
## Data: HH8
##
## REML criterion at convergence: 222.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1357 -0.6398 -0.1682 0.4261 2.4440
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 0.00 0.00
## Residual 62.57 7.91
## Number of obs: 32, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 19.77938 23.68747 30.00000 0.835 0.410
## osmolality_mmol_kg_mean 0.04658 0.06734 30.00000 0.692 0.494
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.998
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 24.554 24.554 1 18.124 0.3924 0.5388HD_CEWL_osml <- lmerTest::lmer(data = HD8,
CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
(1|trial_number))
summary(HD_CEWL_osml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
## Data: HD8
##
## REML criterion at convergence: 171.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55684 -0.59339 -0.07624 0.40035 2.96626
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 6.99 2.644
## Residual 11.28 3.358
## Number of obs: 31, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.11489 11.95116 25.69290 0.177 0.861
## osmolality_mmol_kg_mean 0.06101 0.03291 26.24328 1.854 0.075 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.994## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 31.917 31.917 1 26.042 2.8298 0.1045CH_CEWL_osml <- lmerTest::lmer(data = CH8,
CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
(1|trial_number))
summary(CH_CEWL_osml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
## Data: CH8
##
## REML criterion at convergence: 192.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2869 -0.7878 -0.1247 0.6553 2.0089
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 6.866 2.620
## Residual 30.995 5.567
## Number of obs: 30, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 22.02423 13.70141 27.96594 1.607 0.119
## osmolality_mmol_kg_mean 0.02489 0.03901 27.97274 0.638 0.529
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.993## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 11.376 11.376 1 27.97 0.367 0.5495CD_CEWL_osml <- lmerTest::lmer(data = CD8,
CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
(1|trial_number))
summary(CD_CEWL_osml)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
## Data: CD8
##
## REML criterion at convergence: 170
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8562 -0.5019 0.1381 0.4084 1.8415
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 21.14 4.598
## Residual 11.34 3.368
## Number of obs: 30, groups: trial_number, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.45589 11.16955 27.92145 -0.130 0.8972
## osmolality_mmol_kg_mean 0.07357 0.03133 27.93722 2.348 0.0262 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## osmllty_m__ -0.981## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 55.393 55.393 1 27.937 4.8837 0.03548 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Even though only data for one treatment group was significant, all have a positive relationship
CEWL ~ WVP
The Rsq is better… 0.22 vs 0.15
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_WVP_lm)
## W = 0.98635, p-value = 0.209##
## Call:
## lm(formula = delta_CEWL ~ e_a_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.6389 -6.4501 -0.5578 4.6940 22.9179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7748 1.3846 0.560 0.577
## e_a_mean_tmttrial 3.0979 0.5049 6.136 9.37e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.388 on 131 degrees of freedom
## Multiple R-squared: 0.2232, Adjusted R-squared: 0.2173
## F-statistic: 37.65 on 1 and 131 DF, p-value: 9.373e-09CEWL_WVP_poly <- lm(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ poly(e_a_mean_tmttrial, 2))
plot(CEWL_WVP_poly)
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_WVP_poly)
## W = 0.98651, p-value = 0.2166##
## Call:
## lm(formula = delta_CEWL ~ poly(e_a_mean_tmttrial, 2), data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.3251 -6.0525 -0.4557 5.0184 22.4236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.0035 0.7271 11.008 < 2e-16 ***
## poly(e_a_mean_tmttrial, 2)1 51.4627 8.3851 6.137 9.44e-09 ***
## poly(e_a_mean_tmttrial, 2)2 8.7057 8.3851 1.038 0.301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.385 on 130 degrees of freedom
## Multiple R-squared: 0.2296, Adjusted R-squared: 0.2178
## F-statistic: 19.37 on 2 and 130 DF, p-value: 4.327e-08## Analysis of Variance Table
##
## Model 1: delta_CEWL ~ poly(e_a_mean_tmttrial, 2)
## Model 2: delta_CEWL ~ e_a_mean_tmttrial
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 130 9140.3
## 2 131 9216.1 -1 -75.788 1.0779 0.3011CEWL ~ VPD
We also want to assess whether the vapor pressure deficit lizards experienced during acclimation has an effect on the change in CEWL.
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_VPD_lm)
## W = 0.96688, p-value = 0.002494##
## Call:
## lm(formula = delta_CEWL ~ VPD_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.8041 -5.6085 -0.7356 5.5386 27.3599
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.9609 1.4186 9.841 < 2e-16 ***
## VPD_mean_tmttrial -2.9512 0.5943 -4.965 2.09e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.731 on 131 degrees of freedom
## Multiple R-squared: 0.1584, Adjusted R-squared: 0.152
## F-statistic: 24.66 on 1 and 131 DF, p-value: 2.095e-06Even though the data is slightly nonlinear, a linear model does a fine job explaining the data.
To be extra sure, also check for high leverage and influential points using H-values and Cook’s distance:
# compute values for observations
high_leverage <- data.frame(H = hatvalues(CEWL_VPD_lm)) %>%
mutate(row = row_number())
# compute cutoff value
h_bar <- (3*sum(high_leverage$H))/nrow(high_leverage)
# add to original dataframe
# see which observations have extremely high leverage (if any)
high_leverage_dat <- dat_no_rehab_deltaCEWL %>%
mutate(row = row_number()) %>%
left_join(., high_leverage, by = "row") %>%
dplyr::filter(H > h_bar)
high_leverage_dat## # A tibble: 0 × 38
## # Groups: individual_ID [0]
## # … with 38 variables: measurement_date <date>, type <fct>,
## # individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## # trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## # tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## # CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## # cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## # msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …# get Cook's distance
cooks <- data.frame(c = cooks.distance(CEWL_VPD_lm)) %>%
mutate(row = row_number())
# add to original dataframe
influential <- dat_no_rehab_deltaCEWL %>%
mutate(row = row_number()) %>%
left_join(., cooks, by = "row")
# see moderately influential points
cook_mod_inf <- influential %>%
dplyr::filter(c>0.5)
cook_mod_inf## # A tibble: 0 × 38
## # Groups: individual_ID [0]
## # … with 38 variables: measurement_date <date>, type <fct>,
## # individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## # trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## # tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## # CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## # cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## # msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …NO high leverage or influential points! :)
We will double check a comparison of a polynomial model, just to be sure:
CEWL_VPD_poly <- lm(data = dat_no_rehab_deltaCEWL,
delta_CEWL ~ poly(VPD_mean_tmttrial, 2))
plot(CEWL_VPD_poly)
##
## Shapiro-Wilk normality test
##
## data: residuals(CEWL_VPD_poly)
## W = 0.97275, p-value = 0.008935##
## Call:
## lm(formula = delta_CEWL ~ poly(VPD_mean_tmttrial, 2), data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.2213 -5.2906 -0.9053 4.9843 26.5733
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.0035 0.7523 10.639 < 2e-16 ***
## poly(VPD_mean_tmttrial, 2)1 -43.3514 8.6760 -4.997 1.84e-06 ***
## poly(VPD_mean_tmttrial, 2)2 -14.1269 8.6760 -1.628 0.106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.676 on 130 degrees of freedom
## Multiple R-squared: 0.1752, Adjusted R-squared: 0.1625
## F-statistic: 13.81 on 2 and 130 DF, p-value: 3.647e-06LINE assumptions are equally-well-met.
The polynomial factor is not significant, but the R-sq is slightly higher for the poly model.
Compare WVP / VPD / poly
Compare RMSE and AIC:
## [1] 8.664654## [1] 8.577628## [1] 8.324289## [1] 8.289991CEWL_VPD_models <- list(CEWL_WVP_lm, CEWL_WVP_poly,
CEWL_VPD_lm, CEWL_VPD_poly)
CEWL_VPD_mod_names <- data.frame(model =
c('WVP linear',
'WVP polynomial',
'VPD linear',
'VPD polynomial'
))
CEWL_VPD_AICc <- data.frame(aictab(cand.set = CEWL_VPD_models,
modnames = CEWL_VPD_mod_names))
CEWL_VPD_AICc## model K AICc Delta_AICc ModelLik AICcWt LL
## 1 WVP linear 3 947.3249 0.000000 1.000000000 0.621392996 -470.5695
## 2 WVP polynomial 4 948.3532 1.028205 0.598037054 0.371616037 -470.0203
## 4 VPD polynomial 4 957.4260 10.101091 0.006405839 0.003980544 -474.5568
## 3 VPD linear 3 957.9847 10.659766 0.004844638 0.003010424 -475.8993
## Cum.Wt
## 1 0.6213930
## 2 0.9930090
## 4 0.9969896
## 3 1.0000000RMSE is slightly lower for the polynomial model. But, AIC is not meaningfully different between the two versions for each variable.
I’ll use the lm for each as the best models.
##
## Call:
## lm(formula = delta_CEWL ~ VPD_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.8041 -5.6085 -0.7356 5.5386 27.3599
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.9609 1.4186 9.841 < 2e-16 ***
## VPD_mean_tmttrial -2.9512 0.5943 -4.965 2.09e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.731 on 131 degrees of freedom
## Multiple R-squared: 0.1584, Adjusted R-squared: 0.152
## F-statistic: 24.66 on 1 and 131 DF, p-value: 2.095e-06anova(CEWL_VPD_lm) # single factor slr = can't specify ddf and type SS the same way (& don't need to I think)## Analysis of Variance Table
##
## Response: delta_CEWL
## Df Sum Sq Mean Sq F value Pr(>F)
## VPD_mean_tmttrial 1 1879.3 1879.34 24.656 2.095e-06 ***
## Residuals 131 9985.1 76.22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## lm(formula = delta_CEWL ~ e_a_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.6389 -6.4501 -0.5578 4.6940 22.9179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7748 1.3846 0.560 0.577
## e_a_mean_tmttrial 3.0979 0.5049 6.136 9.37e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.388 on 131 degrees of freedom
## Multiple R-squared: 0.2232, Adjusted R-squared: 0.2173
## F-statistic: 37.65 on 1 and 131 DF, p-value: 9.373e-09## Analysis of Variance Table
##
## Response: delta_CEWL
## Df Sum Sq Mean Sq F value Pr(>F)
## e_a_mean_tmttrial 1 2648.4 2648.41 37.645 9.373e-09 ***
## Residuals 131 9216.1 70.35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1stats for paper
VPD: (estimate = -2.95, SE = 0.59, SS = 1879, F(1,131) = 24.66, p < 0.0001, adj-Rsq = 0.15)
WVP: (estimate = 3.10, SE = 0.50, SS = 2648, F(1,131) = 37.65, p < 0.0001, adj-Rsq = 0.22)
Recovery Models
not in publication
I want to know how the 2-day recovery period affects physiology relative to post- and pre- experiment. To do this, I’ll use two-sided t-tests comparing recovery delta values to the hypothesis of mu=0.
SMI
SMI_rmod_post_exp <- t.test(recovery_v_post_exp$delta_SMI_10_8,
mu = 0, alternative = "two.sided")
SMI_rmod_post_exp##
## One Sample t-test
##
## data: recovery_v_post_exp$delta_SMI_10_8
## t = 6.677, df = 131, p-value = 6.292e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3171205 0.5841444
## sample estimates:
## mean of x
## 0.4506324SMI_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_SMI_10_0,
mu = 0, alternative = "two.sided")
SMI_rmod_pre_exp##
## One Sample t-test
##
## data: recovery_v_pre_exp$delta_SMI_10_0
## t = -11.542, df = 131, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.3426322 -0.9497445
## sample estimates:
## mean of x
## -1.146188Hematocrit
hct_rmod_post_exp <- t.test(recovery_v_post_exp$delta_hct_10_8,
mu = 0, alternative = "two.sided")
hct_rmod_post_exp##
## One Sample t-test
##
## data: recovery_v_post_exp$delta_hct_10_8
## t = -4.2083, df = 126, p-value = 4.85e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -3.033119 -1.092866
## sample estimates:
## mean of x
## -2.062992hct_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_hct_10_0,
mu = 0, alternative = "two.sided")
hct_rmod_pre_exp##
## One Sample t-test
##
## data: recovery_v_pre_exp$delta_hct_10_0
## t = -22.249, df = 129, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -13.50271 -11.29729
## sample estimates:
## mean of x
## -12.4Osmolality
osml_rmod_post_exp <- t.test(recovery_v_post_exp$delta_osml_10_8,
mu = 0, alternative = "two.sided")
osml_rmod_post_exp##
## One Sample t-test
##
## data: recovery_v_post_exp$delta_osml_10_8
## t = 3.4782, df = 121, p-value = 0.0007021
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 4.330203 15.772256
## sample estimates:
## mean of x
## 10.05123osml_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_osml_10_0,
mu = 0, alternative = "two.sided")
osml_rmod_pre_exp##
## One Sample t-test
##
## data: recovery_v_pre_exp$delta_osml_10_0
## t = 3.75, df = 130, p-value = 0.0002651
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 5.758213 18.618378
## sample estimates:
## mean of x
## 12.1883Group Export
Now, I put all the recovery t-test statistics into one dataframe, and export it.
recovery_stats <- broom.mixed::tidy(osml_rmod_pre_exp) %>%
rbind(broom.mixed::tidy(osml_rmod_post_exp)) %>%
rbind(broom.mixed::tidy(hct_rmod_pre_exp)) %>%
rbind(broom.mixed::tidy(hct_rmod_post_exp)) %>%
rbind(broom.mixed::tidy(SMI_rmod_pre_exp)) %>%
rbind(broom.mixed::tidy(SMI_rmod_post_exp)) %>%
mutate(response = c(rep("Plasma Osmolality (mmol/kg)", 2),
rep("Hematocrit (%)", 2),
rep("Body Condition (g')", 2)),
pre_post_exp = c(rep(c("pre", "post"), 3)))
write.csv(recovery_stats,
"./results_statistics/recovery_stats.csv")Figures
Colors & Shapes
Create custom colors, shapes, and labels to make it easy to apply changes across all figures.
CH_color <- brewer.pal(4, "Spectral")[c(4)]
HH_color <- brewer.pal(4, "Spectral")[c(2)]
CD_color <- brewer.pal(4, "Spectral")[c(3)]
HD_color <- brewer.pal(4, "Spectral")[c(1)]
my_colors <- c(CH_color, HH_color, CD_color, HD_color)
CH_shp <- 15
HH_shp <- 19
CD_shp <- 22
HD_shp <- 21
CH_shp_box <- 22
HH_shp_box <- 21
my_shapes <- c(CH_shp, HH_shp, CD_shp, HD_shp)
my_shapes_box <- c(CH_shp_box, HH_shp_box, CD_shp, HD_shp)
my_labels <- c("Cool Humid\n0.6 kPa",
"Hot Humid\n1.1 kPa",
"Cool Dry\n2.5 kPa",
"Hot Dry\n3.8 kPa")Body Mass
Means Only F4
ggplot() +
#plot these first so they end up on the "bottom"
geom_smooth(data = dat_no_rehab,
aes(x = day_n,
y = mass_g,
color = tmt,
group = tmt),
method = "lm",
se = F,
size = 0.7) +
geom_errorbar(data = means,
aes(x = day_n,
y = mean_mass,
color = tmt,
group = tmt,
ymin = mean_mass-se_mass,
ymax = mean_mass+se_mass),
width = .1,
alpha = 0.8) +
geom_point(data = means,
aes(x = day_n,
y = mean_mass,
color = tmt,
#fill = tmt,
shape = tmt),
fill = "white",
alpha = 1,
size = 3) +
theme_classic() +
scale_shape_manual(values = my_shapes, name = "",
labels = my_labels) +
scale_fill_manual(values = my_colors, name = "",
labels = my_labels) +
scale_color_manual(values = my_colors, name = "",
labels = my_labels) +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_y_continuous(breaks = c(seq(9,12)),
labels = c(seq(9,12))) +
xlab("Day") +
ylab("Body Mass (g)") +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
#geom_vline(xintercept = 9,
# linetype = "dashed",
# color = "darkgrey") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 10.8, unit = "pt")
) -> mass_fig_min## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.## `geom_smooth()` using formula = 'y ~ x'
Ending Values F5
ggplot() +
geom_jitter(data = end_vals,
aes(x = tmt,
y = mass_g,
color = tmt,
fill = tmt,
shape = tmt),
size = 1,
alpha = 0.4,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = mass_emmeans,
aes(x = tmt,
y = emmean,
color = tmt,
group = tmt,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 1) +
geom_point(data = mass_emmeans,
aes(x = tmt,
y = emmean,
#color = tmt,
shape = tmt,
fill = tmt),
color = "black",
size = 4) +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "") +
scale_fill_manual(values = my_colors, name = "") +
scale_color_manual(values = my_colors, name = "") +
scale_y_continuous(limits = c(7,15),
breaks = c(seq(7,15, by = 2)),
labels = c(seq(7,15, by = 2))) +
scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
"Hot Humid\n1.1 kPa",
"Cool Dry\n2.5 kPa",
"Hot Dry\n3.8 kPa")) +
xlab("") +
annotate(geom = "text", x = 4, y = 14.4, label = "B",
size = 3) +
annotate(geom = "text", x = 2, y = 15, label = "BC",
size = 3) +
annotate(geom = "text", x = 3, y = 14.6, label = "AC",
size = 3) +
annotate(geom = "text", x = 1, y = 14.8, label = "A",
size = 3) +
ylab("Body Mass (g)") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0, #top
0, #right
0, #bottom
3.4), "mm")
) -> mass_end_boxplot
mass_end_boxplot
## tmt emmean SE df lower.CL upper.CL response
## 1 Cool Humid 10.667647 0.2517552 46.72727 10.161103 11.17419 Body Mass (g)
## 2 Hot Humid 9.727273 0.2559428 47.91893 9.212643 10.24190 Body Mass (g)
## 3 Cool Dry 10.581818 0.2559914 47.91893 10.067090 11.09655 Body Mass (g)
## 4 Hot Dry 9.547059 0.2521465 46.16198 9.039562 10.05456 Body Mass (g)## contrast estimate SE df t.ratio p.value
## 1 Cool Humid - Hot Humid 0.94037433 0.3575715 126.2004 2.6298915 0.04674770
## 2 Cool Humid - Cool Dry 0.08582888 0.3585976 127.4556 0.2393459 0.99515381
## 3 Cool Humid - Hot Dry 1.12058824 0.3554223 126.8612 3.1528361 0.01074007
## 4 Hot Humid - Cool Dry -0.85454545 0.3615808 127.8215 -2.3633597 0.08951418
## 5 Hot Humid - Hot Dry 0.18021390 0.3583741 127.2163 0.5028653 0.95826222
## 6 Cool Dry - Hot Dry 1.03475936 0.3577985 126.4621 2.8920171 0.02307102
## response
## 1 Body Mass (g)
## 2 Body Mass (g)
## 3 Body Mass (g)
## 4 Body Mass (g)
## 5 Body Mass (g)
## 6 Body Mass (g)Hct
Ind + Means
ggplot() +
geom_line(data = dat[complete.cases(dat$hematocrit_percent),],
aes(x = day_n,
y = hematocrit_percent,
color = tmt,
group = individual_ID),
alpha = 0.2) +
geom_line(data = means[complete.cases(means$mean_hct),],
aes(x = day_n,
y = mean_hct,
color = tmt,
group = tmt),
alpha = 1,
size = 1) +
geom_point(data = means,
aes(x = day_n,
y = mean_hct,
color = tmt,
shape = tmt),
alpha = 1,
size = 5) +
geom_vline(xintercept = 9,
linetype = "dashed",
color = "darkgrey") +
theme_classic() +
scale_shape_manual(values = c(15:18), name = "") +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
scale_color_brewer(palette = "Set2", name = "") +
xlab("") +
ylab("Hematocrit (%)") +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 22),
axis.text = element_text(color = "black",
family = "sans",
size = 16),
legend.text = element_text(color = "black",
family = "sans",
size = 22),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0.1, #top
0.1, #right
0.1, #bottom
0.41 #left
), "cm")
) -> hct_fig
hct_fig## Warning: Removed 12 rows containing missing values (`geom_point()`).
Means Only F5
ggplot() +
geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$hematocrit_percent),],
aes(x = day_n,
y = hematocrit_percent,
color = tmt,
group = tmt),
method = "lm",
se = F,
size = 0.7) +
geom_errorbar(data = means[complete.cases(means$mean_hct),],
aes(x = day_n,
y = mean_hct,
color = tmt,
group = tmt,
ymin = mean_hct-se_hct,
ymax = mean_hct+se_hct),
width = .1,
alpha = 0.8) +
geom_point(data = means[complete.cases(means$mean_hct),],
aes(x = day_n,
y = mean_hct,
color = tmt,
#fill = tmt,
shape = tmt),
alpha = 1,
fill = "white",
size = 3) +
theme_classic() +
scale_shape_manual(values = my_shapes, name = "",
labels = my_labels) +
scale_fill_manual(values = my_colors, name = "",
labels = my_labels) +
scale_color_manual(values = my_colors, name = "",
labels = my_labels) +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_y_continuous(breaks = c(25, 30, 35, 40),
labels = c(25, 30, 35, 40),) +
xlab("") +
ylab("Hematocrit (%)") +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
#geom_vline(xintercept = 9,
# linetype = "dashed",
# color = "darkgrey") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 0, l = 10.8, unit = "pt")
) -> hct_fig_min
hct_fig_min## `geom_smooth()` using formula = 'y ~ x'
Ending Values F5
ggplot() +
geom_jitter(data = end_vals,
aes(x = tmt,
y = hematocrit_percent,
color = tmt,
fill = tmt,
shape = tmt),
size = 1,
alpha = 0.4,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = hct_emmeans,
aes(x = tmt,
y = emmean,
color = tmt,
group = tmt,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 1) +
geom_point(data = hct_emmeans,
aes(x = tmt,
y = emmean,
#color = tmt,
shape = tmt,
fill = tmt),
color = "black",
size = 4) +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "") +
scale_fill_manual(values = my_colors, name = "") +
scale_color_manual(values = my_colors, name = "") +
xlab("") +
scale_y_continuous(limits = c(10,55),
breaks = c(seq(10,50, by = 10)),
labels = c(seq(10,50, by = 10))) +
scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
"Hot Humid\n1.1 kPa",
"Cool Dry\n2.5 kPa",
"Hot Dry\n3.8 kPa")) +
annotate(geom = "text", x = 4, y = 47.6, label = "AB",
size = 3) +
annotate(geom = "text", x = 2, y = 44.6, label = "B",
size = 3) +
annotate(geom = "text", x = 3, y = 54.6, label = "A",
size = 3) +
annotate(geom = "text", x = 1, y = 54.6, label = "A",
size = 3) +
ylab("Hematocrit (%)") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0, #top
0, #right
0, #bottom
2.24), "mm")
) -> hct_end_boxplot
hct_end_boxplot## Warning: Removed 3 rows containing missing values (`geom_point()`).
## tmt emmean SE df lower.CL upper.CL response
## 1 Cool Humid 30.33352 1.808392 5.860156 25.88283 34.78422 Hematocrit (%)
## 2 Hot Humid 25.77326 1.830146 6.142261 21.32008 30.22645 Hematocrit (%)
## 3 Cool Dry 29.70335 1.815476 5.950456 25.25206 34.15464 Hematocrit (%)
## 4 Hot Dry 27.63092 1.815279 5.948653 23.17978 32.08205 Hematocrit (%)## contrast estimate SE df t.ratio p.value
## 1 Cool Humid - Hot Humid 4.5602630 1.272328 123.0398 3.5841894 0.002703439
## 2 Cool Humid - Cool Dry 0.6301769 1.254531 123.1123 0.5023206 0.958383983
## 3 Cool Humid - Hot Dry 2.7026074 1.253006 123.0665 2.1568990 0.141276348
## 4 Hot Humid - Cool Dry -3.9300861 1.288296 123.2298 -3.0506087 0.014665591
## 5 Hot Humid - Hot Dry -1.8576556 1.283706 123.1026 -1.4471033 0.472669199
## 6 Cool Dry - Hot Dry 2.0724305 1.262089 123.0607 1.6420637 0.359065321
## response
## 1 Hematocrit (%)
## 2 Hematocrit (%)
## 3 Hematocrit (%)
## 4 Hematocrit (%)
## 5 Hematocrit (%)
## 6 Hematocrit (%)Osml
Ind + Means
ggplot() +
geom_line(data = dat[complete.cases(dat$osmolality_mmol_kg_mean),],
aes(x = day_n,
y = osmolality_mmol_kg_mean,
color = tmt,
group = individual_ID),
alpha = 0.2) +
geom_line(data = means[complete.cases(means$mean_osml),],
aes(x = day_n,
y = mean_osml,
color = tmt,
group = tmt),
alpha = 1,
size = 1) +
geom_point(data = means,
aes(x = day_n,
y = mean_osml,
color = tmt,
shape = tmt),
alpha = 1,
size = 5) +
geom_vline(xintercept = 9,
linetype = "dashed",
color = "darkgrey") +
scale_shape_manual(values = c(15:18), name = "") +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
ylim(300,450) +
scale_color_brewer(palette = "Set2", name = "") +
xlab("") +
ylab("Plasma Osmolality (mmol/kg)") +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 22),
axis.text = element_text(color = "black",
family = "sans",
size = 16),
legend.text = element_text(color = "black",
family = "sans",
size = 22),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0.6, #top
0.1, #right
0.1, #bottom
0.1 #left
), "cm")
) -> osml_fig
osml_fig## Warning: Removed 1 row containing missing values (`geom_line()`).## Warning: Removed 12 rows containing missing values (`geom_point()`).
Means Only F4
ggplot() +
geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$osmolality_mmol_kg_mean),],
aes(x = day_n,
y = osmolality_mmol_kg_mean,
color = tmt,
group = tmt),
method = "lm",
se = F,
size = 0.7) +
geom_errorbar(data = means,
aes(x = day_n,
y = mean_osml,
color = tmt,
group = tmt,
ymin = mean_osml-se_osml,
ymax = mean_osml+se_osml),
width = .1,
alpha = 0.8) +
geom_point(data = means,
aes(x = day_n,
y = mean_osml,
color = tmt,
#fill = tmt,
shape = tmt),
fill = "white",
alpha = 1,
size = 3) +
theme_classic() +
scale_shape_manual(values = my_shapes, name = "",
labels = my_labels) +
scale_fill_manual(values = my_colors, name = "",
labels = my_labels) +
scale_color_manual(values = my_colors, name = "",
labels = my_labels) +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_y_continuous(breaks = c(seq(340,380, by = 10)),
labels = c(seq(340,380, by = 10))) +
xlab("") +
ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
#geom_vline(xintercept = 9,
# linetype = "dashed",
# color = "darkgrey") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 0, l = 1, unit = "pt")
) -> osml_fig_min
osml_fig_min## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 12 rows containing missing values (`geom_point()`).
Stats! Check Pairwise Diffs ~ Time
Since Plasma osmolality has such a nonlinear trend, we need to test whether the elevated values in the middle of the experiment are significantly different than the values taken before and/or after.
# first make sub-dfs for each tmt group
HH <- dat_no_rehab %>%
dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD <- dat_no_rehab %>%
dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH <- dat_no_rehab %>%
dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD <- dat_no_rehab %>%
dplyr::filter(substr(tmt, 1, 6) == "Cool D")
# next do pairwise tests for osml on the diff exp days, for each tmt group
pair_HH <- TukeyHSD(aov(data = HH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_HD <- TukeyHSD(aov(data = HD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CH <- TukeyHSD(aov(data = CH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CD <- TukeyHSD(aov(data = CD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
# put into a df and export
osml_pairwise_df <- as.data.frame(pair_HD[[1]]) %>%
rbind(as.data.frame(pair_HH[[1]])) %>%
rbind(as.data.frame(pair_CD[[1]])) %>%
rbind(as.data.frame(pair_CH[[1]])) %>%
mutate(day_diff = paste("day", substr(rownames(.), 1, 3)),
tmt = c(rep("Hot Dry",3),
rep("Hot Humid",3),
rep("Cool Dry",3),
rep("Cool Humid",3)),
CI_95 = paste(round(lwr, digits = 2), round(upr, digits = 2), sep = ", "),
diff = round(diff, digits = 2)) %>%
dplyr::select(tmt, day_diff, diff, CI_95, p_adj = "p adj")
# save
write.csv(osml_pairwise_df, "./results_statistics/osmolality_pairwise_diffs.csv")Nope, none of the differences between days within tmt groups are significantly different.
Ending Values F5
ggplot() +
geom_jitter(data = end_vals,
aes(x = tmt,
y = osmolality_mmol_kg_mean,
color = tmt,
fill = tmt,
shape = tmt),
size = 1,
alpha = 0.4,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = osml_emmeans,
aes(x = tmt,
y = emmean,
color = tmt,
group = tmt,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 1) +
geom_point(data = osml_emmeans,
aes(x = tmt,
y = emmean,
#color = tmt,
shape = tmt,
fill = tmt),
color = "black",
size = 4) +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "") +
scale_fill_manual(values = my_colors, name = "") +
scale_color_manual(values = my_colors, name = "") +
scale_y_continuous(limits = c(290,470),
breaks = c(seq(300,450, by = 50)),
labels = c(seq(300,450, by = 50))) +
scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
"Hot Humid\n1.1 kPa",
"Cool Dry\n2.5 kPa",
"Hot Dry\n3.8 kPa")) +
xlab("") +
annotate(geom = "text", x = 4, y = 427, label = "A",
size = 3) +
annotate(geom = "text", x = 2, y = 452, label = "A",
size = 3) +
annotate(geom = "text", x = 3, y = 437, label = "A",
size = 3) +
annotate(geom = "text", x = 1, y = 470, label = "A",
size = 3) +
ylab(bquote('Osmolality (mmol '*kg^-1*')')) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0, #top
0, #right
0, #bottom
0), "mm")
) -> osml_end_boxplot
osml_end_boxplot## Warning: Removed 10 rows containing missing values (`geom_point()`).
## tmt emmean SE df lower.CL upper.CL
## 1 Cool Humid 350.3202 10.16423 4.874923 323.9893 376.6511
## 2 Hot Humid 354.0963 10.13690 4.824373 327.7506 380.4420
## 3 Cool Dry 349.1611 10.17830 4.903510 322.8413 375.4808
## 4 Hot Dry 362.0244 10.15727 4.863220 335.6918 388.3569
## response
## 1 Plasma Osmolality (mmol/kg)
## 2 Plasma Osmolality (mmol/kg)
## 3 Plasma Osmolality (mmol/kg)
## 4 Plasma Osmolality (mmol/kg)## contrast estimate SE df t.ratio p.value
## 1 Cool Humid - Hot Humid -3.776107 5.031118 116.0115 -0.7505502 0.87626200
## 2 Cool Humid - Cool Dry 1.159115 5.128594 116.0589 0.2260104 0.99590771
## 3 Cool Humid - Hot Dry -11.704175 5.078569 116.0347 -2.3046207 0.10288296
## 4 Hot Humid - Cool Dry 4.935222 5.085161 116.0412 0.9705144 0.76645671
## 5 Hot Humid - Hot Dry -7.928068 5.035558 116.0197 -1.5744172 0.39720576
## 6 Cool Dry - Hot Dry -12.863291 5.114607 116.0115 -2.5150105 0.06282759
## response
## 1 Plasma Osmolality (mmol/kg)
## 2 Plasma Osmolality (mmol/kg)
## 3 Plasma Osmolality (mmol/kg)
## 4 Plasma Osmolality (mmol/kg)
## 5 Plasma Osmolality (mmol/kg)
## 6 Plasma Osmolality (mmol/kg)CEWL
Ind + Means
ggplot() +
geom_line(data = dat[complete.cases(dat$CEWL_g_m2h_mean),],
aes(x = day_n,
y = CEWL_g_m2h_mean,
color = tmt,
group = individual_ID),
alpha = 0.2) +
geom_line(data = means[complete.cases(means$mean_CEWL),],
aes(x = day_n,
y = mean_CEWL,
color = tmt,
group = tmt),
alpha = 1,
size = 1) +
geom_point(data = means[complete.cases(means$mean_CEWL),],
aes(x = day_n,
y = mean_CEWL,
color = tmt,
shape = tmt),
alpha = 1,
size = 5) +
theme_classic() +
scale_shape_manual(values = c(15:18), name = "") +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_color_brewer(palette = "Set2", name = "") +
xlab("Day") +
ylim(0,60) +
ylab(bquote('CEWL (g/'*m^2*'h)')) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 22),
axis.text = element_text(color = "black",
family = "sans",
size = 16),
legend.text = element_text(color = "black",
family = "sans",
size = 22),
legend.text.align = 0,
legend.position = "bottom"
) -> CEWL_fig
CEWL_fig
Means Only F2
ggplot() +
geom_errorbar(data = means[complete.cases(means$mean_CEWL),],
aes(x = day_n,
y = mean_CEWL,
color = tmt,
group = tmt,
ymin = mean_CEWL-se_CEWL,
ymax = mean_CEWL+se_CEWL),
width = .1,
#position=position_dodge(.1),
alpha = 0.8) +
geom_line(data = means[complete.cases(means$mean_CEWL),],
aes(x = day_n,
y = mean_CEWL,
color = tmt,
#linetype = tmt,
group = tmt),
alpha = 0.7,
size = .5) +
geom_point(data = means[complete.cases(means$mean_CEWL),],
aes(x = day_n,
y = mean_CEWL,
color = tmt,
shape = tmt),
fill = "white",
alpha = 1,
size = 3) +
theme_classic() +
scale_shape_manual(values = my_shapes, name = "",
labels = my_labels) +
scale_fill_manual(values = my_colors, name = "",
labels = my_labels) +
scale_color_manual(values = my_colors, name = "",
labels = my_labels) +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_y_continuous(breaks = c(20, 25, 30, 35, 40),
limits = c(18,40)) +
xlab("Day") +
ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none"
#legend.position = c(0.25,0.85)
) -> CEWL_fig_min
CEWL_fig_min
# use ggarrange so legend is centered
CEWL_fig_formatted <- ggarrange(CEWL_fig_min,
ncol = 1, nrow = 1,
common.legend = TRUE,
legend = "bottom")
# save
ggsave(filename = "experiment_CEWL_fig.pdf",
plot = CEWL_fig_formatted,
path = "./results_figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 90)delta ~ VPD F3
ggplot(data = dat_no_rehab_deltaCEWL) +
geom_point(aes(x = VPD_mean_tmttrial,
y = delta_CEWL,
color = tmt,
fill = tmt,
shape = tmt),
size = 2,
alpha = 0.4
) +
geom_smooth(aes(x = VPD_mean_tmttrial,
y = delta_CEWL),
se = F,
formula = y ~ x,
method = "lm",
color = "black") +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_fill_manual(values = my_colors, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_color_manual(values = my_colors, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_y_continuous(limits = c(-13,40),
breaks = c(seq(-10,40, by = 10)),
labels = c(seq(-10,40, by = 10))
) +
scale_x_continuous(limits = c(0,5),
breaks = c(0.6, 1.1, 2.5, 3.8),
labels = c("0.6\nCH",
"1.1\nHH",
"2.5\nCD",
"3.8\nHD")) +
xlab("Vapor Pressure Deficit (kPa)") +
ylab(expression(Delta ~ 'CEWL')) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "bottom",
legend.spacing.y = unit(0, "mm")
#plot.margin = unit(c(0, #top
# 0, #right
# 0, #bottom
# 0), "mm")
) -> CEWL_VPD_fig
CEWL_VPD_fig
# code to export this figure alone
# NA bc now grouped with the next figure
# # use ggarrange so legend is centered
# CEWL_VPD_fig_formatted <- ggarrange(CEWL_VPD_fig,
# ncol = 1, nrow = 1,
# common.legend = TRUE,
# legend = "bottom")
# # save
# ggsave(filename = "exp_CEWL_delta_VPD.pdf",
# plot = CEWL_VPD_fig_formatted,
# path = "./results_figures",
# device = "pdf",
# dpi = 600,
# units = "mm",
# width = 80, height = 90)delta ~ WVP F3
ggplot(data = dat_no_rehab_deltaCEWL) +
geom_point(aes(x = e_a_mean_tmttrial,
y = delta_CEWL,
color = tmt,
fill = tmt,
shape = tmt),
size = 2,
alpha = 0.4
) +
geom_smooth(aes(x = e_a_mean_tmttrial,
y = delta_CEWL),
se = F,
formula = y ~ x,
method = "lm",
color = "black") +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_fill_manual(values = my_colors, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_color_manual(values = my_colors, name = "",
labels = c("Cool Humid (CH)",
"Hot Humid (HH)",
"Cool Dry (CD)",
"Hot Dry (HD)")) +
scale_y_continuous(limits = c(-13,40),
breaks = c(seq(-10,40, by = 10)),
labels = c(seq(-10,40, by = 10))
) +
scale_x_continuous(limits = c(0,5),
breaks = c(0.5, 2.0, 2.3, 4.5),
labels = c("0.5\nCD",
"2.0 \nHD ",
" 2.3\n CH",
"4.5\nHH")) +
xlab("Water Vapor Pressure (kPa)") +
ylab(expression(Delta ~ 'CEWL')) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "bottom",
legend.spacing.y = unit(0, "mm")
#plot.margin = unit(c(0, #top
# 0, #right
# 0, #bottom
# 0), "mm")
) -> CEWL_WVP_fig
CEWL_WVP_fig
Ending Values F5
ggplot() +
geom_jitter(data = end_vals,
aes(x = tmt,
y = CEWL_g_m2h_mean,
color = tmt,
fill = tmt,
shape = tmt),
size = 1,
alpha = 0.4,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = CEWL_emmeans,
aes(x = tmt,
y = emmean,
color = tmt,
group = tmt,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 1) +
geom_point(data = CEWL_emmeans,
aes(x = tmt,
y = emmean,
#color = tmt,
shape = tmt,
fill = tmt),
color = "black",
size = 4) +
theme_classic() +
scale_shape_manual(values = my_shapes_box, name = "") +
scale_fill_manual(values = my_colors, name = "") +
scale_color_manual(values = my_colors, name = "") +
scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
"Hot Humid\n1.1 kPa",
"Cool Dry\n2.5 kPa",
"Hot Dry\n3.8 kPa")) +
scale_y_continuous(limits = c(9,63),
breaks = c(seq(10,60, by = 10)),
labels = c(seq(10,60, by = 10))) +
xlab("") +
ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
annotate(geom = "text", x = 4, y = 50, label = "C",
size = 3) +
annotate(geom = "text", x = 2, y = 63, label = "B",
size = 3) +
annotate(geom = "text", x = 3, y = 41, label = "C",
size = 3) +
annotate(geom = "text", x = 1, y = 52, label = "A",
size = 3) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text.align = 0,
legend.position = "none",
plot.margin = unit(c(0, #top
0, #right
0, #bottom
0), "mm")
) -> CEWL_end_boxplot
CEWL_end_boxplot## Warning: Removed 1 rows containing missing values (`geom_point()`).
## tmt emmean SE df lower.CL upper.CL response
## 1 Cool Humid 30.36157 1.445712 10.68500 27.16810 33.55505 CEWL (g/m2h)
## 2 Hot Humid 36.76428 1.446643 10.69577 33.56916 39.95941 CEWL (g/m2h)
## 3 Cool Dry 23.72439 1.446548 10.69551 20.52946 26.91931 CEWL (g/m2h)
## 4 Hot Dry 24.61233 1.434910 10.37491 21.43074 27.79392 CEWL (g/m2h)## contrast estimate SE df t.ratio p.value
## 1 Cool Humid - Hot Humid -6.4027092 1.477747 125.0660 -4.3327496 1.740823e-04
## 2 Cool Humid - Cool Dry 6.6371867 1.482776 125.4155 4.4761900 9.849699e-05
## 3 Cool Humid - Hot Dry 5.7492451 1.468760 125.2015 3.9143543 8.436015e-04
## 4 Hot Humid - Cool Dry 13.0398959 1.482912 125.4316 8.7934385 9.026113e-14
## 5 Hot Humid - Hot Dry 12.1519543 1.469693 125.2717 8.2683648 1.057265e-12
## 6 Cool Dry - Hot Dry -0.8879416 1.467098 125.0811 -0.6052366 9.302438e-01
## response
## 1 CEWL (g/m2h)
## 2 CEWL (g/m2h)
## 3 CEWL (g/m2h)
## 4 CEWL (g/m2h)
## 5 CEWL (g/m2h)
## 6 CEWL (g/m2h)Exp CEWL ~ Osml
end_vals_CEWL_osml <- dat %>%
dplyr::filter(day_n == 8) %>%
dplyr::filter(complete.cases(CEWL_g_m2h_mean, osmolality_mmol_kg_mean))
ggplot(end_vals_CEWL_osml) +
aes(x = osmolality_mmol_kg_mean,
y = CEWL_g_m2h_mean,
color = tmt,
shape = tmt) +
geom_point(size = 1.5,
alpha = 0.8) +
stat_smooth(formula = y ~ x,
method = "lm",
se = F,
size = 1,
alpha = 0.9) +
theme_classic() +
xlab(bquote('Osmolality (mmol '*kg^-1*')')) +
ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) +
scale_shape_manual(values = c(21,19, 22,15), name = "") +
scale_fill_brewer(palette = "Spectral", name = "") +
scale_color_brewer(palette = "Spectral", name = "") +
scale_x_continuous(breaks = c(300, 350, 400, 450)) +
scale_y_continuous(breaks = c(20, 30, 40, 50),
limits = c(12,57)) +
guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.position = "bottom"
) -> exp_end_CEWL_osml_fig
exp_end_CEWL_osml_fig
Multi-Figures
hydration ~ time F4
ggarrange(osml_fig_min,
hct_fig_min,
mass_fig_min,
ncol = 1, nrow = 3,
labels = c("A", "B", "C"),
common.legend = TRUE,
legend = "bottom"
) -> experiment_multi_fig## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 12 rows containing missing values (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 12 rows containing missing values (`geom_point()`).## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'end values F5
ggarrange(CEWL_end_boxplot,
osml_end_boxplot,
hct_end_boxplot,
mass_end_boxplot,
ncol = 2, nrow = 2,
labels = c("A", "B", "C", "D"),
widths = c(2, 2.045), heights = c(2, 2),
common.legend = FALSE
) -> ending_values_multi_fig## Warning: Removed 1 rows containing missing values (`geom_point()`).## Warning: Removed 10 rows containing missing values (`geom_point()`).## Warning: Removed 3 rows containing missing values (`geom_point()`).VPD/WVP F3
# use ggarrange so legend is centered
CEWL_VPD_fig_formatted <- ggarrange(CEWL_VPD_fig,
CEWL_WVP_fig,
ncol = 1, nrow = 2,
common.legend = TRUE,
legend = "bottom",
labels = c("A", "B"))
# save
ggsave(filename = "exp_CEWL_delta_VPD_WVP.pdf",
plot = CEWL_VPD_fig_formatted,
path = "./results_figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 150)